{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Binned Likelihood Tutorial\n",
    "\n",
    "The detection, flux determination, and spectral modeling of Fermi LAT sources is accomplished by a maximum likelihood optimization technique as described in the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/) (see also, e.g., [Abdo, A. A. et al. 2009, ApJS, 183, 46](http://adsabs.harvard.edu/abs/2009ApJS..183...46A)).\n",
    "\n",
    "To illustrate how to use the Likelihood software, this tutorial gives a step-by-step description for performing a binned likelihood analysis."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Binned vs Unbinned Likelihood\n",
    "\n",
    "Binned likelihood analysis is the preferred method for most types of LAT analysis (see [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/)).\n",
    "\n",
    "However, when analyzing data over short time periods (with few events), it is better to use the **unbinned** analysis.\n",
    "\n",
    "To perform an unbinned likelihood analysis, see the [Unbinned Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html) tutorial.\n",
    "\n",
    "**Additional references**:\n",
    "* [SciTools References](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/references.html)\n",
    "* Descriptions of available [Spectral and Spatial Models](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html)\n",
    "* Examples of [XML Model Definitions for Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#xmlModelDefinitions):\n",
    "    * [Power Law](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#powerlaw)\n",
    "    * [Broken Power Law](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#brokenPowerLaw)\n",
    "    * [Broken Power Law 2](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#powerLaw2)\n",
    "    * [Log Parabola](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#logParabola)\n",
    "    * [Exponential Cutoff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#expCutoff)\n",
    "    * [BPL Exponential Cutoff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#bplExpCutoff)\n",
    "    * [Gaussian](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#gaussian)\n",
    "    * [Constant Value](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#constantValue)\n",
    "    * [File Function](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#fileFunction)\n",
    "    * [Band Function](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#bandFunction)\n",
    "    * [PL Super Exponential Cutoff](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html#plSuperExpCutoff)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Prerequisites\n",
    "\n",
    "You will need an **event** data file, a **spacecraft** data file (also referred to as the \"pointing and livetime history\" file), and the current **background models** (available for [download](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html)). They are also found in code cells below.\n",
    "\n",
    "You may choose to select your own data files, or to use the files provided within this tutorial.\n",
    "\n",
    "Custom data sets may be retrieved from the [Lat Data Server](http://fermi.gsfc.nasa.gov/cgi-bin/ssc/LAT/LATDataQuery.cgi)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Outline\n",
    "\n",
    "1. **Make Subselections from the Event Data**\n",
    "\n",
    "    Since there is computational overhead for each event associated with each diffuse component, it is useful to filter out any events that are not within the extraction region used for the analysis.\n",
    "\n",
    "\n",
    "2. **Make Counts Maps from the Event Files**\n",
    "    \n",
    "    By making simple FITS images, we can inspect our data and pick out obvious sources.\n",
    "\n",
    "\n",
    "3. **Download the latest diffuse models**\n",
    "\n",
    "    The recommended models for a normal point source analysis are `gll_iem_v07.fits` (a very large file) and `iso_P8R3_SOURCE_V2_v1.txt`. All of the background models along with a description of the models are available [here](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html).\n",
    "\n",
    "\n",
    "4. **Create a Source Model XML File**\n",
    "\n",
    "    The source model XML file contains the various sources and their model parameters to be fit using the **gtlike** tool.\n",
    "\n",
    "\n",
    "5. **Create a 3D Counts Cube**\n",
    "\n",
    "    The binned counts cube is used to reduce computation requirements in regions with large numbers of events.\n",
    "\n",
    "\n",
    "6. **Compute Livetimes**\n",
    "\n",
    "    Precomputing the livetime for the dataset speeds up the exposure calculation.\n",
    "    \n",
    "\n",
    "7. **Compute Exposure Cube**\n",
    "\n",
    "    This accounts for exposure as a function of energy, based on the cuts made. The exposure map must be recomputed if any change is made to the data selection or binning.\n",
    "    \n",
    "\n",
    "8. **Compute Source Maps**\n",
    "\n",
    "    Here the exposure calculation is applied to each of the sources described in the model.\n",
    "\n",
    "\n",
    "9. **Perform the Likelihood Fit**\n",
    "\n",
    "    Fitting the data to the model provides flux, errors, spectral indices, and other information.\n",
    "\n",
    "\n",
    "10. **Create a Model Map**\n",
    "\n",
    "    This can be compared to the counts map to verify the quality of the fit and to make a residual map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Make subselections from the event data\n",
    "\n",
    "For this case we will use two years of LAT Pass 8 data. This is a longer data set than is described in the [Extract LAT Data](../DataSelection/1.ExtractLATData.ipynb) tutorial.\n",
    "\n",
    ">**NOTE**: The ROI used by the binned likelihood analysis is defined by the 3D counts map boundary. The region selection used in the data extraction step, which is conical, must fully contain the 3D counts map spatial boundary, which is square.\n",
    "\n",
    "Selection of data:\n",
    "\n",
    "    Search Center (RA, DEC) =(193.98, -5.82)\n",
    "    Radius = 15 degrees\n",
    "    Start Time (MET) = 239557417 seconds (2008-08-04 T15:43:37)\n",
    "    Stop Time (MET) = 302572802 seconds (2010-08-04 T00:00:00)\n",
    "    Minimum Energy = 100 MeV\n",
    "    Maximum Energy = 500000 MeV\n",
    "\n",
    "This two-year dataset generates numerous data files. We provide the user with the original event data files and the accompanying spacecraft file:\n",
    "\n",
    "* L181126210218F4F0ED2738_PH00.fits (5.0 MB)\n",
    "* L181126210218F4F0ED2738_PH01.fits (10.5 MB)\n",
    "* L181126210218F4F0ED2738_PH02.fits (6.5 MB)\n",
    "* L181126210218F4F0ED2738_PH03.fits (9.2 MB)\n",
    "* L181126210218F4F0ED2738_PH04.fits (7.4 MB)\n",
    "* L181126210218F4F0ED2738_PH05.fits (6.2 MB)\n",
    "* L181126210218F4F0ED2738_PH06.fits (4.5 MB)\n",
    "* L181126210218F4F0ED2738_SC00.fits (256 MB spacecraft file)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH00.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH01.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH02.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH03.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH04.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH05.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_PH06.fits\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/L181126210218F4F0ED2738_SC00.fits"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mkdir ./data\n",
    "!mv *.fits ./data\n",
    "!ls ./data"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to combine the two events files for your analysis, you must first generate a text file listing the events files to be included.\n",
    "\n",
    "If you do not wish to download all the individual files, you can skip to the next step and retrieve the combined, filtered event file. However, you will need the spacecraft file to complete the analysis, so you should retrieve that now.\n",
    "\n",
    "To generate the file list, type:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!ls ./data/*_PH* > ./data/binned_events.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When analyzing point sources, it is recommended that you include events with high probability of being photons. To do this, you should use **gtselect** to cut on the event class, keeping only the SOURCE class events (event class 128, or as recommended in the Cicerone).\n",
    "\n",
    "In addition, since we do not wish to cut on any of the three event types (conversion type, PSF, or EDISP), we will use `evtype=3` (which corresponds to standard analysis in Pass 7). Note that `INDEF` is the default for evtype in gtselect."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "gtselect evclass=128 evtype=3\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Be aware that `evclass` and `evtype` are hidden parameters. So, to use them, you must type them on the command line.\n",
    "\n",
    "The text file you made (`binned_events.txt`) will be used in place of the input fits filename when running gtselect. The syntax requires that you use an @ before the filename to indicate that this is a text file input rather than a fits file.\n",
    "\n",
    "We perform a selection to the data we want to analyze. For this example, we consider the source class photons within our 15 degree region of interest (ROI) centered on the blazar 3C 279. For some of the selections that we made with the data server and don't want to modify, we can use \"INDEF\" to instruct the tool to read those values from the data file header. Here, we are only filtering on event class (not on event type) and applying a zenith cut, so many of the parameters are designated as \"INDEF\"."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We apply the **gtselect** tool to the data file as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtselect evclass=128 evtype=3\n",
    "    @./data/binned_events.txt\n",
    "    ./data/3C279_binned_filtered.fits\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    INDEF\n",
    "    100\n",
    "    500000\n",
    "    90"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the last step we also selected the energy range and the maximum zenith angle value (90 degrees) as suggested in Cicerone and recommended by the LAT instrument team.\n",
    "\n",
    "The Earth's limb is a strong source of background gamma rays and we can filter them out with a zenith-angle cut. The use of \"zmax\" in calculating the exposure allows for a more selective method than just using the ROI cuts in controlling the Earth limb contamination. The filtered data from the above steps are provided [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_filtered.fits)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After the data selection is made, we need to select the good time intervals in which the satellite was working in standard data taking mode and the data quality was good. For this task we use **gtmktime** to select GTIs by filtering on information provided in the spacecraft file. The current **gtmktime** filter expression recommended by the LAT team in the Cicerone is:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "(DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This excludes time periods when some spacecraft event has affected the quality of the data; it ensures the LAT instrument was in normal science data-taking mode.\n",
    "\n",
    "Here is an example of running **gtmktime** for our analysis of the region surrounding 3C 279."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtmktime\n",
    "    @./data/L181126210218F4F0ED2738_SC00.fits\n",
    "    (DATA_QUAL>0)&&(LAT_CONFIG==1)\n",
    "    no\n",
    "    ./data/3C279_binned_filtered.fits\n",
    "    ./data/3C279_binned_gti.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The data file with all the cuts described above is provided in this [link](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_gti.fits). A more detailed discussion of data selection can be found in the [Data Preparation](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data_preparation.html) analysis thread.\n",
    "\n",
    "To view the DSS keywords in a given extension of a data file, use the **gtvcut** tool and review the data cuts on the EVENTS extension. This provides a listing of the keywords reflecting each cut applied to the data file and their values, including the entire list of GTIs. (Use the option `suppress_gtis=no` to view the entire list.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtvcut suppress_gtis=no\n",
    "    ./data/3C279_binned_gti.fits\n",
    "    EVENTS"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here you can see the event class and event type, the location and radius of the data selection, as well as the energy range in MeV, the zenith angle cut, and the fact that the time cuts to be used in the exposure calculation are defined by the GTI table.\n",
    "\n",
    "Various Fermitools will be unable to run if you have multiple copies of a particular DSS keyword. This can happen if the position used in extracting the data from the data server is different than the position used with **gtselect**. It is wise to review the keywords for duplicates before proceeding. If you do have keyword duplication, it is advisable to regenerate the data file with consistent cuts."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Make a counts map from the event data\n",
    "\n",
    "Next, we create a counts map of the ROI, summed over photon energies, in order to identify candidate sources and to ensure that the field looks sensible as a simple sanity check. For creating the counts map, we will use the [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) tool with the option \"CMAP\" (no spacecraft file is necessary for this step).\n",
    "\n",
    "Then we will view the output file, as shown below:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/3C279_binned_gti.fits\n",
    "    ./data/3C279_binned_cmap.fits\n",
    "    NONE\n",
    "    150\n",
    "    150\n",
    "    0.2\n",
    "    CEL\n",
    "    193.98\n",
    "    -5.82\n",
    "    0.0\n",
    "    AIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We chose an ROI of 15 degrees, corresponding to 30 degrees in diameter. Since we want a pixel size of 0.2 degrees/pixel, then we must select 30/0.2=150 pixels for the size of the x and y axes. The last command launches the visualization tool _ds9_ and produces a display of the generated [counts](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_cmap.fits) map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/BinnedLikelihood/3C279_binned_counts_map.png'>\n",
    "\n",
    "You can see several strong sources and a number of weaker sources in this map. Mousing over the positions of these sources shows that two of them are likely 3C 279 and 3C 273.\n",
    "\n",
    "It is important to inspect your data prior to proceeding to verify that the contents are as you expect. A malformed data query or improper data selection can generate a non-circular region, or a file with zero events. By inspecting your data prior to analysis, you have an opportunity to detect such issues early in the analysis.\n",
    "\n",
    "A more detailed discussion of data exploration can be found in the [Explore LAT Data](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/explore_latdata.html) analysis thread."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Create a 3-D (binned) counts map\n",
    "\n",
    "Since the counts map shows the expected data, you are ready to prepare your data set for analysis. For binned likelihood analysis, the data input is a three-dimensional counts map with an energy axis, called a counts cube. The gtbin tool performs this task as well by using the `CCUBE` option.\n",
    "\n",
    "<img src=\"https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/BinnedLikelihood/square_in_circle.png\">\n",
    "\n",
    "The binning of the counts map determines the binning of the exposure calculation. The likelihood analysis may lose accuracy if the energy bins are not sufficiently narrow to accommodate more rapid variations in the effective area with decreasing energy below a few hundred MeV. For a typical analysis, ten logarithmically spaced bins per decade in energy are recommended. The analysis is less sensitive to the spatial binning and 0.2 deg bins are a reasonable standard.\n",
    "\n",
    "This counts cube is a square binned region that must fit within the circular acceptance cone defined during the data extraction step, and visible in the counts map above. To find the maximum size of the region your data will support, find the side of a square that can be fully inscribed within your circular acceptance region (multiply the radius of the acceptance cone by sqrt(2)). For this example, the maximum length for a side is 21.21 degrees."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To create the counts cube, we run [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) as follows:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CCUBE\n",
    "    ./data/3C279_binned_gti.fits\n",
    "    ./data/3C279_binned_ccube.fits\n",
    "    NONE\n",
    "    100\n",
    "    100\n",
    "    0.2\n",
    "    CEL\n",
    "    193.98\n",
    "    -5.82\n",
    "    0.0\n",
    "    AIT\n",
    "    LOG\n",
    "    100\n",
    "    500000\n",
    "    37"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "[gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) takes the following as parameters:\n",
    "* Type of output file (CCUBE|CMAP|LC|PHA1|PHA2|HEALPIX)\n",
    "* Event data file name\n",
    "* Output file name\n",
    "* Spacecraft data file name\n",
    "* Size of the X axis in pixels\n",
    "* Size of the Y axis in pixels\n",
    "* Image scale (in degrees/pixel)\n",
    "* Coordindate system (CEL - celestial; GAL - galactic) (pick CEL or GAL)\n",
    "* First coordinate of image center in degrees (RA or galactic l)\n",
    "* Second coordinate of image center in degrees (DEC or galactic b)\n",
    "* Rotation angle of image axis, in degrees\n",
    "* Projection method (AIT|ARC|CAR|GLS|MER|NCP|SIN|STG|TAN)\n",
    "* Algorithm for defining energy bins (FILE|LIN|LOG)\n",
    "* Start value for first energy bin in MeV\n",
    "* Stop value for last energy bin in MeV\n",
    "* Number of logarithmically uniform energy bins"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The counts cube generated in this step is provided [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_ccube.fits).\n",
    "\n",
    "If you open the file with _ds9_, you see that it is made up of 37 images, one for each logarithmic energy bin. By playing through these images, it is easy to see how the PSF of the LAT changes with energy. You can also see that changing energy cuts could be helpful when trying to optimize the localization or spectral information for specific sources.\n",
    "\n",
    "Be sure to verify that there are no black corners on your counts cube. These corners correspond to regions with no data and will cause errors in your exposure calculations."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 4. Download the latest diffuse model files\n",
    "\n",
    "When you use the current Galactic diffuse emission model ([`gll_iem_v07.fits`](https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/gll_iem_v07.fits)) in a likelihood analysis, you also want to use the corresponding model for the extragalactic isotropic diffuse emission, which includes the residual cosmic-ray background. The recommended isotropic model for point source analysis is [`iso_P8R3_SOURCE_V2_v1.txt`](https://fermi.gsfc.nasa.gov/ssc/data/analysis/software/aux/4fgl/iso_P8R3_SOURCE_V2_v1.txt).\n",
    "\n",
    "All the Pass 8 background models have been included in the Fermitools distribution, in the `$(FERMI_DIR)/refdata/fermi/galdiffuse/` directory. If you use that path in your model, you should not have to download the diffuse models individually.\n",
    "\n",
    ">**NOTE**: Keep in mind that the isotropic model needs to agree with both the event class and event type selections you are using in your analysis. The iso_P8R3_SOURCE_V2_v1.txt isotropic spectrum is valid only for the latest response functions and only for data sets with front + back events combined. All of the most up-to-date background models along with a description of the models are available [here](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/BackgroundModels.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 5. Create a source model XML file\n",
    "\n",
    "The [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) tool reads the source model from an XML file. The model file contains your best guess at the locations and spectral forms for the sources in your data. A source model can be created using the [model editor](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/modeleditor.txt) tool, by using the user contributed tool `make4FGLxml.py` (available at the [user-contributed tools](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/) page), or by editing the file directly within a text editor.\n",
    "\n",
    "Here we cannot use the same source model that was used to analyze six months of data in the Unbinned Likelihood tutorial, as the 2-year data set contains many more significant sources and will not converge. Instead, we will use the 4FGL catalog to define our source model by running `make4FGLxml.py`. To run the script, you will need to download the current LAT catalog file and place it in your working directory:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/make4FGLxml.py\n",
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/gll_psc_v18.fit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!mv make4FGLxml.py gll_psc_v18.fit ./data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!python ./data/make4FGLxml.py ./data/gll_psc_v18.fit ./data/3C279_binned_gti.fits -o ./data/3C279_input_model.xml "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Note that we are using a high level of significance so that we only fit the brightest sources, and we have forced the extended sources to be modeled as point sources.\n",
    "\n",
    "It is also necessary to specify the entire path to location of the diffuse model on your system. Clearly, the simple 4-source model we used for the 6-month [Unbinned Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html) analysis would have been too simplistic."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This XML file uses the spectral model from the 4FGL catalog analysis for each source. (The catalog file is available at the [LAT 8-yr Catalog page](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/).) However, that analysis used a subset of the available spectral models. A dedicated analysis of the region may indicate a different spectral model is preferred.\n",
    "\n",
    "For more details on the options available for your XML models, see:\n",
    "* Descriptions of available [Spectral and Spatial Models](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/source_models.html)\n",
    "* Examples of [XML Model Definitions for Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/xml_model_defs.html)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, the `make4FGLxml.py` script automatically adds 10 degrees to your ROI to account for sources that lie outside your data region, but which may contribute photons to your data. In addition, it gives you the ability to free only some of the spectral parameters for sources within your ROI, and fixes them for the others.\n",
    "\n",
    "With hundreds of sources, there are too many free parameters to gain a good spectral fit. It is advisable to revise these values so that only sources near your source of interest, or very bright source, have all spectral parameters free. Farther away, you can fix the spectral form and free only the normalization parameter (or \"prefactor\"). If you are working in a crowded region or have nested sources (e.g. a point source on top of an extended source), you will probably want to fix parameters for some sources even if they lie close to your source of interest.\n",
    "\n",
    "Only the normalization parameter will be left free for the remaining sources within the ROI. We have also used the significance parameter (`-s`) of `make4FLGxml.py` to free only the brightest sources in our ROI. In addition, we used the `-v` flag to override that for sources that are significantly variable. Both these changes are necessary: having too many free parameters will not allow the fit to converge (see the section for the fitting step)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### XML for Extended Sources\n",
    "\n",
    "In some regions, the [make4FGLxml.py](https://fermi.gsfc.nasa.gov/ssc/data/analysis/user/make4FGLxml.py) script may add one or more extended sources to your XML model.\n",
    "\n",
    "The script will provide the number of extended sources included in the model. In order to use these extended sources, you will need to downloaded the extended source templates from the [LAT Catalog](https://fermi.gsfc.nasa.gov/ssc/data/access/lat/8yr_catalog/) page (look for \"Extended Source template archive\").\n",
    "\n",
    "Extract the archive in the directory of your choice and note the path to the template files, which have names like `W44.fits` and `VelaX.fits`. You will need to provide the path to the template file to the script before you run it.\n",
    "\n",
    "Here is an example of the proper format for an extended source XML entry for Binned Likelihood analysis:\n",
    "\n",
    "```xml\n",
    "<source name=\"SpatialMap_source\" type=\"DiffuseSource\">\n",
    "<spectrum type=\"PowerLaw2\">\n",
    "<parameter free=\"1\" max=\"1000.0\" min=\"1e-05\" name=\"Integral\" scale=\"1e-06\" value=\"1.0\"/>\n",
    "<parameter free=\"1\" max=\"-1.0\" min=\"-5.0\" name=\"Index\" scale=\"1.0\" value=\"-2.0\"/>\n",
    "<parameter free=\"0\" max=\"200000.0\" min=\"20.0\" name=\"LowerLimit\" scale=\"1.0\" value=\"20.0\"/>\n",
    "<parameter free=\"0\" max=\"200000.0\" min=\"20.0\" name=\"UpperLimit\" scale=\"1.0\" value=\"2e5\"/>\n",
    "</spectrum>\n",
    "<spatialModel W44 file=\"$(PATH_TO_FILE)/W44.fits\" type=\"SpatialMap\" map_based_integral=\"true\">\n",
    "<parameter free=\"0\" max=\"1000.0\" min=\"0.001\" name=\"Normalization\" scale= \"1.0\" value=\"1.0\"/>\n",
    "</spatialModel>\n",
    "</source>\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 6. Compute livetimes and exposure\n",
    "\n",
    "To speed up the exposure calculations performed by Likelihood, it is helpful to pre-compute the livetime as a function of sky position and off-axis angle. The [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gtltcube.txt) tool creates a livetime cube, which is a [HealPix](http://healpix.jpl.nasa.gov/) table, covering the entire sky, of the integrated livetime as a function of inclination with respect to the LAT z-axis.\n",
    "\n",
    "Here is an example of how to run [gtltcube](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/gtltcube.txt):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtltcube zmax=90\n",
    "    ./data/3C279_binned_gti.fits\n",
    "    ./data/L181126210218F4F0ED2738_SC00.fits\n",
    "    ./data/3C279_binned_ltcube.fits\n",
    "    0.025\n",
    "    1"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    ">**Note**: Values such as \"0.1\" for \"Step size in cos(theta) are known to give unexpected results. Use \"0.09\" instead."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The livetime cube generated from this analysis can be found [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_ltcube.fits).\n",
    "\n",
    "For more information about the livetime cubes see the documentation in the [Cicerone](https://fermi.gsfc.nasa.gov/ssc/data/analysis/documentation/Cicerone/Cicerone_Likelihood/) and also the explanation in the [Unbinned Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html) tutorial."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 7. Compute exposure map\n",
    "\n",
    "Next, you must apply the livetime calculated in the previous step to your region of interest. To do this, we use the [gtexpcube2](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpcube2.txt) tool, which is an updated version of the previous **gtexpcube**. This tool generates a binned exposure map, an accounting of the exposure at each position in the sky, that are a required input to the likelihood process.\n",
    "\n",
    ">**NOTE**: In the past, running **gtsrcmaps** calculated the exposure map for you, so most analyses skipped the binned exposure map generation step. With the introduction of **gtexpcube2**, this is no longer the case. You must explicitly command the creation of the exposure map as a separate analysis step."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to create an exposure map that accounts for contributions from all the sources in your analysis region, you must consider not just the sources included in the counts cube. The large PSF of the LAT means that at low energies, sources from well outside your counts cube could affect the sources you are analyzing. To compensate for this, you must create an exposure map that includes sources up to 10 degrees outside your ROI. (The ROI is determined by the radius you downloaded from the data server, here a 15 degree radius.) In addition, you should account for all the exposure that contributes to those additional sources. Since the exposure map uses square pixels, to match the binning in the counts cube, and to ensure we don't have errors, we generate a 300x300 pixel map."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you provide [gtexpcube2](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtexpcube2.txt) a filename for your counts cube, it will use the information from that file to define the geometry of the exposure map. This is legacy behavior and will not give you the necessary 20° buffer you need to completely account for the exposure of nearby sources. (It will also cause an error in the next step.)\n",
    "\n",
    "Instead, you should specify the appropriate geometry for the exposure map, remembering that the counts cube used 0.2 degree pixel binning. To do that, enter `none` when asked for a Counts cube.\n",
    "\n",
    "**Note**: If you get a \"`File not found`\" error in the examples below, just put the IRF name in explicitly. The appropriate IRF for this data set is `P8R3_SOURCE_V2`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtexpcube2\n",
    "    ./data/3C279_binned_ltcube.fits\n",
    "    none\n",
    "    ./data/3C279_binned_expcube.fits\n",
    "    P8R3_SOURCE_V2\n",
    "    300\n",
    "    300\n",
    "    .2\n",
    "    193.98\n",
    "    -5.82\n",
    "    0\n",
    "    AIT\n",
    "    CEL\n",
    "    100\n",
    "    500000\n",
    "    37"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The generated exposure map can be found [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_expcube.fits).\n",
    "\n",
    "At this point, you may decide it is easier to simply generate exposure maps for the entire sky. You may be right, as it certainly simplifies the step when scripting. However, making an all-sky map increases the processing time for this step, though the increase is modest.\n",
    "\n",
    "To generate an all-sky exposure map (rather than the exposure map we calculated above) you need to specify the proper binning and explicitly give the number of pixels for the entire sky (360°x180°).\n",
    "\n",
    "Here is an example:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtexpcube2\n",
    "    ./data/3C279_binned_ltcube.fits\n",
    "    none\n",
    "    ./data/3C279_binned_allsky_expcube.fits\n",
    "    P8R3_SOURCE_V2\n",
    "    1800\n",
    "    900\n",
    "    .2\n",
    "    193.98\n",
    "    -5.82\n",
    "    0\n",
    "    AIT\n",
    "    CEL\n",
    "    100\n",
    "    500000\n",
    "    37"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The all-sky exposure map can be found [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_allsky_expcube.fits).\n",
    "\n",
    "Just as in the [Unbinned Likelihood](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/likelihood_tutorial.html) analysis, the exposure needs to be recalculated if the ROI, zenith angle, time, event class, or energy selections applied to the data are changed. For the binned analysis, this also includes the spatial and energy binning of the 3D counts map (which affects the exposure map as well)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 8. Compute source map\n",
    "\n",
    "The [gtsrcmaps](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtsrcmaps.txt) tool creates model counts maps for use with the binned likelihood analysis. To do this, it takes each source spectrum in the XML model, multiplies it by the exposure at the source position, and convolves that exposure with the effective PSF.\n",
    "\n",
    "This is an example of how to run the tool:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtsrcmaps\n",
    "    ./data/3C279_binned_ltcube.fits\n",
    "    ./data/3C279_binned_ccube.fits\n",
    "    ./data/3C279_input_model.xml\n",
    "    ./data/3C279_binned_allsky_expcube.fits\n",
    "    ./data/3C279_binned_srcmaps.fits\n",
    "    CALDB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The output file from [gtsrcmaps](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtsrcmaps.txt) can be found [here](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_srcmaps.fits).\n",
    "\n",
    "Because your model map can include sources outside your ROI, you may see a list of warnings at the beginning of the output. These are expected (because you have properly included sources outside your ROI in your XML file) and should cause no problem in your analysis. In addition, if your exposure map is too small for the region, you will see the following warning:\n",
    "\n",
    "```\n",
    "Caught St13runtime_error at the top level:\n",
    "Request for exposure at a sky position that is outside of the map boundaries.\n",
    "\n",
    "The contribution of the diffuse source outside of the exposure\n",
    "and counts map boundaries is being computed to account for PSF\n",
    "leakage into the analysis region. To handle this, use an all-sky\n",
    "binned exposure map. Alternatively, to neglect contributions\n",
    "outside of the counts map region, use the emapbnds=no option when\n",
    "running gtsrcmaps.\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this situation, you should increase the dimensions of your exposure map, or just move to the all-sky version.\n",
    "\n",
    "Source map generation for the point sources is fairly quick, and maps for many point sources may take up a lot of disk space. If you are analyzing a single long data set, it may be preferable to pre-compute only the source maps for the diffuse components at this stage.\n",
    "\n",
    "[gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) will compute maps for the point sources on the fly if they appear in the XML definition and a corresponding map is not in the source maps FITS file. To skip generating source maps for point sources, specify \"`ptsrc=no`\" on the command line when running **gtsrcmaps**. However, if you expect to perform multiple fits on the same set of data, precomputing the source maps will probably save you time."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9. Run gtlike\n",
    "\n",
    ">NOTE: Prior to running **gtlike** for Unbinned Likelihood, it is necessary to calculate the diffuse response for each event (when that response is not precomputed). However, for Binned Likelihood analysis the diffuse response is calculated over the entire bin, so this step is not necessary.\n",
    "\n",
    "If you want to use the **energy dispersion correction** during your analysis, you must enable this feature using the environment variable `USE_BL_EDISP`. This may be set on the command line using:\n",
    "\n",
    "```bash\n",
    "export USE_BL_EDISP=true\n",
    "```\n",
    "or, depending on your shell,\n",
    "\n",
    "```\n",
    "setenv USE_BL_EDISP=true\n",
    "```\n",
    "\n",
    "To disable the use of energy dispersion, you must unset the variable:\n",
    "```bash\n",
    "unset USE_BL_EDISP\n",
    "```\n",
    "or\n",
    "```\n",
    "unsetenv USE_BL_EDISP\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```bash\n",
    "export USE_BL_EDISP=true\n",
    "```\n",
    "or, depending on your shell,\n",
    "\n",
    "```\n",
    "setenv USE_BL_EDISP=true\n",
    "```\n",
    "\n",
    "To disable the use of energy dispersion, you must unset the variable:\n",
    "```bash\n",
    "unset USE_BL_EDISP\n",
    "```\n",
    "or\n",
    "```\n",
    "unsetenv USE_BL_EDISP\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we are ready to run the [gtlike](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtlike.txt) application.\n",
    "\n",
    "Here, we request that the fitted parameters be saved to an output XML model file for use in later steps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!wget https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_output_model.xml"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtlike refit=yes plot=yes sfile=./data/3C279_binned_output.xml\n",
    "    BINNED\n",
    "    ./data/3C279_binned_srcmaps.fits\n",
    "    ./data/3C279_binned_allsky_expcube.fits\n",
    "    ./data/3C279_binned_ltcube.fits\n",
    "    ./data/3C279_input_model.xml\n",
    "    CALDB\n",
    "    NEWMINUIT"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Most of the entries prompted for are fairly obvious. In addition to the various XML and FITS files, the user is prompted for a choice of IRFs, the type of statistic to use, and the optimizer."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The statistics available are:\n",
    "* **UNBINNED**: This should be used for short timescale or low source count data. If this option is chosen then parameters for the spacecraft file, event file, and exposure file must be given. See explanation in: [Likelihood Tutorial]()\n",
    "\n",
    "* **BINNED**: This is a standard binned analysis as described in this tutorial. This analysis is used for long timescale or high-density data (such as in the Galactic plane) which can cause memory errors in the unbinned analysis. If this option is chosen then parameters for the source map file, livetime file, and exposure file must be given."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There are five optimizers from which to choose: `DRMNGB`, `DRMNFB`, `NEWMINUIT`, `MINUIT` and `LBFGS`. Generally speaking, the faster way to find the parameter estimates is to use `DRMNGB` (or `DRMNFB`) to find initial values and then use `MINUIT` (or `NEWMINUIT`) to find more accurate results. If you have trouble achieving convergence at first, you can loosen your tolerance by setting the hidden parameter `ftol` on the command line. (The default value for `ftol` is `0.001`.)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Analyzing a 2-year dataset will take many hours (in our case more than 2 days with a 32-bit machine with 1 GB of RAM). The required running time is high if your source is in the Galactic plane. Here is some output from our fit, where 4FGL J1229.0+0202 and 4FGL J1256.1-0547 corresponds to 3C 273 and 3C 279, respectively:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "This is gtlike version \n",
    "\n",
    "...\n",
    "\n",
    "Photon fluxes are computed for the energy range 100 to 500000 MeV\n",
    "\n",
    "4FGL J1229.0+0202:\n",
    "norm: 8.16706 +/- 0.0894921\n",
    "alpha: 2.49616 +/- 0.015028\n",
    "beta: 0.104635 +/- 0.0105201\n",
    "Eb: 279.04\n",
    "TS value: 32017.6\n",
    "Flux: 6.69253e-07 +/- 7.20102e-09 photons/cm^2/s\n",
    "\n",
    "4FGL J1256.1-0547:\n",
    "norm: 2.38177 +/- 0.0296458\n",
    "alpha: 2.25706 +/- 0.0116212\n",
    "beta: 0.0665607 +/- 0.00757385\n",
    "Eb: 442.052\n",
    "TS value: 29261.7\n",
    "Flux: 5.05711e-07 +/- 6.14833e-09 photons/cm^2/s\n",
    "\n",
    "...\n",
    "\n",
    "gll_iem_v07:\n",
    "Prefactor: 0.900951 +/- 0.0235397\n",
    "Index: 0\n",
    "Scale: 100\n",
    "Flux: 0.000469334 +/- 1.22608e-05 photons/cm^2/s\n",
    "\n",
    "iso_P8R3_SOURCE_V2_v1:\n",
    "Normalization: 1.13545 +/- 0.0422581\n",
    "Flux: 0.000139506 +/- 5.19439e-06 photons/cm^2/s\n",
    "\n",
    "WARNING: Fit may be bad in range [100, 199.488] (MeV)\n",
    "WARNING: Fit may be bad in range [251.124, 316.126] (MeV)\n",
    "WARNING: Fit may be bad in range [6302.3, 7933.61] (MeV)\n",
    "WARNING: Fit may be bad in range [39744.4, 50032.1] (MeV)\n",
    "WARNING: Fit may be bad in range [315519, 397190] (MeV)\n",
    "\n",
    "Total number of observed counts: 207751\n",
    "Total number of model events: 207407\n",
    "\n",
    "-log(Likelihood): 73014.38504\n",
    "\n",
    "Writing fitted model to 3C279_binned_output.xml\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Since we selected `plot=yes` in the command line, a plot of the fitted data appears.\n",
    "\n",
    "<img src=\"https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/BinnedLikelihood/3C279_binned_spectral_fit.png\">"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the first plot, the counts/MeV vs MeV are plotted. The points are the data, and the lines are the models. Error bars on the points represent sqrt(Nobs) in that band, where Nobs is the observed number of counts. The black line is the sum of the models for all sources. \n",
    "\n",
    "The colored lines follow the sources as follows:\n",
    "* Black - summed model\n",
    "* Red - first source (see below)\n",
    "* Green - second source\n",
    "* Blue - third source\n",
    "* Magenta - fourth source\n",
    "* Cyan - the fifth source\n",
    "\n",
    "If you have more sources, the colors are reused in the same order. In our case we have, in order of decreasing value on the y-axis: summed model (black), the extragalactic background (black), the galactic background (cyan), 3C 273 (red), and 3C 279 (black).\n",
    "\n",
    "The second plot gives the residuals between your model and the data. Error bars here represent (sqrt(Nopbs))/Npred, where Npred is the predicted number of counts in each band based on the fitted model."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To assess the quality of the fit, look first for the words at the top of the output `<Optimizer> did successfully converge.` Successful convergence is a minimum requirement for a good fit.\n",
    "\n",
    "Next, look at the energy ranges that are generating warnings of bad fits. If any of these ranges affect your source of interest, you may need to revise the source model and refit. You can also look at the residuals on the plot (bottom panel). If the residuals indicate a poor fit overall (e.g., the points trending all low or all high) you should consider changing your model file, perhaps by using a different source model definition, and refit the data.\n",
    "    \n",
    "If the fits and spectral shapes are good, but could be improved, you may wish to simply update your model file to hold some of the spectral parameters fixed. For example, by fixing the spectral model for 3C 273, you may get a better quality fit for 3C 279. Close the plot and you will be asked if you wish to refit the data."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```\n",
    "Refit? [y] n\n",
    "Elapsed CPU time: 1571.805872\n",
    "```"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here, hitting `return` will instruct the application to fit again. We are happy with the result, so we type `n` and end the fit."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Results\n",
    "\n",
    "When it completes, **gtlike** generates a standard output XML file. If you re-run the tool in the same directory, these files will be overwritten by default. Use the `clobber=no` option on the command line to keep from overwriting the output files.\n",
    "\n",
    "Unfortunately, the fit details and the value for the `-log(likelihood)` are not recorded in the automatic output files.\n",
    "\n",
    "You should consider logging the output to a text file for your records by using `> fit_data.txt` (or something similar) with your **gtlike** command.\n",
    "\n",
    "Be aware, however, that this will make it impossible to request a refit when the likelihood process completes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!gtlike plot=yes sfile=./data/3C279_output_model.xml > fit_data.txt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In this example, we used the `sfile` parameter to request that the model results be written to an output XML file. This file contains the source model results that were written to `results.dat` at the completion of the fit.\n",
    "\n",
    ">    **Note**: If you have specified an output XML model file and you wish to modify your model while waiting at the `Refit? [y]` prompt, you will need to copy the results of the output model file to your input model before making those modifications. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The results of the likelihood analysis have to be scaled by the quantity called \"scale\" in the XML model in order to obtain the total photon flux (photons cm-2 s-1) of the source. You must refer to the model formula of your source for the interpretation of each parameter. In our example the 'prefactor' of our power law model of the first fitted source (4FGLJ1159.5-0723) has to be scaled by the factor 'scale'=10-14. For example the total flux of 4FGLJ1159.5-0723 is the integral between 100 MeV and 500000 MeV of:"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "$Prefactor \\cdot scale \\cdot (E /100)^{index}=(6.7017x10-14) \\cdot (E/100)^{-2.0196}$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Errors reported with each value in the `results.dat` file are 1Ïƒ estimates (based on inverse-Hessian at the optimum of the log-likelihood surface)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Other Useful Hidden Parameters\n",
    "\n",
    "If you are scripting and wish to generate multiple output files without overwriting, the `results` and `specfile` parameters allow you to specify output filenames for the `results.dat` and `counts_spectra.fits` files respectively.\n",
    "\n",
    "If you do not specify a source model output file with the `sfile` parameter, then the input model file will be overwritten with the latest fit. This is convenient as it allows the user to edit that file while the application is waiting at the `Refit? [y]` prompt so that parameters can be adjusted and set free or fixed. This would be similar to the use of the \"newpar\", \"freeze\", and \"thaw\" commands in [XSPEC](http://heasarc.gsfc.nasa.gov/docs/xanadu/xspec/index.html). "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 10. Create a model map\n",
    "\n",
    "For comparison to the counts map data, we create a model map of the region based on the fit parameters.\n",
    "\n",
    "This map is essentially an infinite-statistics counts map of the region-of-interest based on our model fit.\n",
    "\n",
    "The [gtmodel](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtmodel.txt) application reads in the fitted model, applies the proper scaling to the source maps, and adds them together to get the final map. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtmodel\n",
    "    ./data/3C279_binned_srcmaps.fits\n",
    "    ./data/3C279_binned_output.xml\n",
    "    ./data/3C279_model_map.fits\n",
    "    CALDB\n",
    "    ./data/3C279_binned_ltcube.fits\n",
    "    ./data/3C279_binned_allsky_expcube.fits"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To understand how well the fit matches the data, we want to compare the [model map](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_model_map.fits) just created with the counts map over the same field of view. First we have to create the [new counts map](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_binned_cmap_small.fits) that matches in size the model map (the one generated in encircles the ROI, while the model map is completely inscribed within the ROI): We will use again the [gtbin](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/help/gtbin.txt) tool with the option `CMAP` as shown below: "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "%%bash\n",
    "gtbin\n",
    "    CMAP\n",
    "    ./data/3C279_binned_gti.fits\n",
    "    ./data/3C279_binned_cmap_small.fits\n",
    "    NONE\n",
    "    100\n",
    "    100\n",
    "    0.2\n",
    "    CEL\n",
    "    193.98\n",
    "    -5.82\n",
    "    0.0\n",
    "    STG"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here we've plotted the model map next to the the energy-summed counts map for the data.\n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/BinnedLikelihood/3C279_binned_map_comparison.png'>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally we want to create the [residual map](https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/data/BinnedLikelihood/3C279_residual.fits) by using the FTOOL **farith** to check if we can improve the model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "%%bash\n",
    "farith\n",
    "    ./data/3C279_binned_cmap_small.fits\n",
    "    ./data/3C279_model_map.fits\n",
    "    ./data/3C279_residual.fits\n",
    "    SUB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The residual map is shown below. As you can see, the binning we chose probably used pixels that were too large.\n",
    "\n",
    "The primary sources, 3C 273 and 3C 279, have some positive pixels next to some negative ones. This effect could be lessened by either using a smaller pixel size or by offsetting the central position slightly from the position of the blazar (or both).\n",
    "\n",
    "If your residual map contains bright sources, the next step would be to iterate the analysis with the additional sources included in the XML model file. \n",
    "\n",
    "<img src='https://fermi.gsfc.nasa.gov/ssc/data/analysis/scitools/images/BinnedLikelihood/3C279_binned_residuals.png'>"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.14"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
